/*
 * KineticsSolver.cpp
 *
 *  Created on: 28 Jul 2011
 *      Author: Allan
 */

#include "KineticsSolver.h"

KineticsSolver::KineticsSolver()
{}

KineticsSolver::KineticsSolver(const Multiphase& multiphase, const ReactionSystem& eReactions, const ReactionSystem& kReactions) :
multiphase(multiphase), eReactions(eReactions), kReactions(kReactions)
{
	// Create a species map object to find the indexes of the species
	const SpeciesMap speciesMap(multiphase.GetSpecies());
	
	// Set the indexes data of the system of equilibrium- and kinetically-controlled reactions
	this->eReactions.SetIndexes(speciesMap);
	this->kReactions.SetIndexes(speciesMap);
	
	// The primary (jSpecies), equilibrium (eSpecies) and kinetic (kSpecies) species
	const vector<string> jSpecies = this->eReactions.Products();
	const vector<string> eSpecies = this->eReactions.Reactants();
	const vector<string> kSpecies = this->kReactions.Reactants();
	
	// Initialize the stoichiometric matrices of the system of equilibrium- and kinetically-controlled reactions
	vej = this->eReactions.StoichiometricMatrix(jSpecies);
	vkj = this->kReactions.StoichiometricMatrix(jSpecies);
	
	// Set the indexes of the primary (jIndexes), equilibrium (eIndexes) and kinetic (kIndexes) species
	jIndexes = speciesMap[jSpecies];
	eIndexes = speciesMap[eSpecies];
	kIndexes = speciesMap[kSpecies];
	
	// Set the number of primary (Nj), equilibrium (Ne) and kinetic (Nk) species
	Nj = jIndexes.size();
	Ne = eIndexes.size();
	Nk = kIndexes.size();
	
	// Initialise the auxiliar Eigen vectors with zero entries
	nj   = VectorXd::Zero(Nj);
	ne   = VectorXd::Zero(Ne);
	nk   = VectorXd::Zero(Nk);
	uj   = VectorXd::Zero(Nj); 
	ujnk = VectorXd::Zero(Nj+Nk);
	
	// Initialise the ODE solver
	odeSolver.Initialiaze(Nj+Nk);
	
	// Initialiaze the equilibrium solver
	eSolver.Initialiaze(multiphase, eReactions);
	eSolver.SetTotalCompositionMode();
	eSolver.UseReducedFormulation();
}

void KineticsSolver::SetControlSolution(double epsRel)
{
	odeSolver.SetControlSolution(0, epsRel);
}

const bool KineticsSolver::Iterate(double& t, double tEnd, double& dt, double T, double P, VectorXd& n)
{
	// Transfer the nj, ne, and nk part of the vector (n) to the vectors (nj), (ne), and (nk) respectively
	for(uint j = 0; j < Nj; ++j) nj[j] = n[jIndexes[j]];
	for(uint e = 0; e < Ne; ++e) ne[e] = n[eIndexes[e]];
	for(uint k = 0; k < Nk; ++k) nk[k] = n[kIndexes[k]];
	
	// Calculate the uj and nk part of the vector (ujnk) = [uj nk], where uj = nj + \sum_e vej * ne
	ujnk.segment(00, Nj) = nj + vej.transpose() * ne;
	ujnk.segment(Nj, Nk) = nk;
	
	// The kinetics problem, which contains the definition of the ODE function to advance the solution in time 
	KineticsProblem problem(*this, T, P, n);
	
	// Use the ODE solver to calculate the vector (ujnk) at the new time (tn) = (t) + (dt)
	odeSolver.Iterate(problem, t, tEnd, dt, ujnk);
	
	// Transfer the uj and nk part of the vector (ujnk) = [uj nk] to the vectors (uj) and (nk) respectively
	uj = ujnk.segment(00, Nj);
	nk = ujnk.segment(Nj, Nk);
	
	// Update the nk part of the vector (n)
	for(uint k = 0; k < Nk; ++k) n[kIndexes[k]] = nk[k];
	
	// Set the total composition of the primary species in the equilibrium solver
	eSolver.SetTotalComposition(uj);
	
	// Solve the equilibrium problem with temperature (T), pressure (P), and composition vector (n), where the nj and ne part of the composition vector (n) will be determined 
	auto diagnostic = eSolver.Solve(T, P, n);
	
	if(!diagnostic.converged)
		cerr << "Error: The equilibrium calculation did not converged." << endl;
	
	return true;
}

const bool KineticsSolver::Solve(double& t, double tEnd, double& dt, double T, double P, VectorXd& n)
{
	while(t < tEnd)
		if(!Iterate(t, tEnd, dt, T, P, n))
			return false;
	
	return true;
}

KineticsProblem::KineticsProblem(KineticsSolver& kSolver, double T, double P, const VectorXd& n) :
kSolver(kSolver), T(T), P(P), n(n), a(n)
{
	Rk = VectorXd::Zero(kSolver.Nk);
}

void KineticsProblem::Function(double t, const VectorXd& ujnk, VectorXd& R)
{
	// Get the necessary references to the data in the kinetic solver (kSolver) to calculate the vector (R)
	    const Multiphase& multiphase = kSolver.multiphase;
	const ReactionSystem& kReactions = kSolver.kReactions;
	          const uint& Nj         = kSolver.Nj;
	          const uint& Nk         = kSolver.Nk;
	 const vector<Index>& kIndexes   = kSolver.kIndexes;
	      const MatrixXd& vkj        = kSolver.vkj;
	   EquilibriumSolver& eSolver    = kSolver.eSolver;
	            VectorXd& uj         = kSolver.uj;
	            VectorXd& nk         = kSolver.nk;
	
	// Transfer the uj and nk part of the vector (ujnk) = [uj nk] to the vectors (uj) and (nk) respectively
	uj = ujnk.segment(00, Nj);
	nk = ujnk.segment(Nj, Nk);
	
	// Update the nk part of the vector (n)
	for(uint k = 0; k < Nk; ++k) n[kIndexes[k]] = nk[k];
	
	// Set the total composition of the primary species in the equilibrium solver
	eSolver.SetTotalComposition(uj);
	
	// Solve the equilibrium problem with temperature (T), pressure (P), and composition vector (n), where the nj and ne part of the composition vector (n) will be determined 
	eSolver.Solve(T, P, n);
	
	// Calculate the activities of the species, now that the composition vector (n) has been determined
	a = multiphase.Activities(T, P, n);
	
	// Calculate the mass rates of the kinetic-controlled reactions with temperature (T), pressure (P), composition vector (n), and activity vector (a)
	Rk = kReactions.MassRates(T, P, n, a);
	
	// Calculate the right-hand side of the ordinary differential equations
	R.segment(00, Nj) = -vkj.transpose() * Rk; // The rates of the total composition variables uj
	R.segment(Nj, Nk) = Rk;                    // The rates of the kinetic-controlled reactions
}
